function [LL, varargout] = loglikelihood(modelId, data, params, takeNegative)
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	%%%%% Inputs:
	% modelId:			       integer (1 => BinomialLogistic, 2 => Poisson)
	% data:			           object:
	% 	.dims:				     	object:
	%		.dims1:						cell(NumParts,1)
	%			{ii}:						integer: gives dim1_ii that corresponds to Xparts{ii}
	%		.Xpart_2_NumX:				1 x NumParts: gives integers (dim2 of Xparts{ii}.X)
	%		.Xpart_2_NumX_FEs:			1 x NumParts: gives integers (dim2 of Xparts{ii}.X_FEs)
	%		.Xpart_2_Num_FE_vals:		cell(1,NumParts)
	%			{ii}:						1 x NumX_FEs_i  --> gives integer: number of possible values for Xparts{ii}.X_FEs(:,ff)
	%		.NumFEvals2Keep:			cell(1,NumParts) --> gives integer
	%		.NumParts
	%		.NumObs
	%		.NumParams
	%    .LL_offset:               1 x 1
	%	 .Xparts:				   cell(NumParts,1)
	%	 	{ii}:				       object
	%			.X:					       dim1_ii x NumX
	%			.X_FEs:					   dim1_ii x NumX_FEs
	%			.NumX_FE_vals:			   1 x NumX_FEs: gives integer
	%    .N:                       NumObs x 1    (if modelId==1)
	%    .logLambdaOffset:         NumObs x 1     (if modelId==2 or modelId==3)
	%    .Y:                       NumObs x 1
	%    .spellDurations:          NumObs x 1 (if modelId==3)
	% params:                  NumX x 1
	% takeNegative:            boolean
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	%%%%% Outputs:
	% LL:            		scalar
	% grad:					NumParams x 1
	% hessian:				NumParams x NumParams
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	
	dims = data.dims;
	
	% Split params into parts
	params_parts = split_combine_parts(dims, 1, params);
	
	% Compute V
	V = compute_Xbeta2(dims, data.Xparts, params_parts);	% NumObs x 1
	
	%%% Compute LL, dLL_dV and d2LL_dV2 (model specific)
	if modelId == 1  % Logistic-Binomial
		p  = 1./(1+exp(-V));								% NumObs x 1
		LL = data.LL_offset + data.Y'*V + data.N'*log(1-p);	% 1 x 1
		clear V;
	end
	if modelId == 2  % Poisson
		logLambda  = data.logLambdaOffset + V;								% NumObs x 1
		clear V;
		lambda     = exp(logLambda);										% NumObs x 1
		LL         = data.LL_offset  + data.Y'*logLambda -sum(lambda,1);	% 1 x 1
		clear logLambda;
	end
	if modelId == 3  % Exponential
		logLambda  = data.logLambdaOffset + V;						% NumObs x 1
		clear V;
		D_lambda = data.spellDurations.*exp(logLambda);				% NumObs x 1
		LL         = data.LL_offset  + data.Y'*logLambda - sum(D_lambda);	% 1 x 1
		clear logLambda;
	end
	
	%%% Output LL and derivatives
	if takeNegative; LL = -LL; end;
	
	if nargout >= 2
		%% Compute dLL_dV
		if modelId == 1  % Logistic-Binomial
			dLL_dV   = data.Y - data.N .* p;	% NumObs x 1
		end
		if modelId == 2  % Poisson
			dLL_dV = data.Y - lambda; 			% NumObs x 1
		end
		if modelId == 3  % Exponential
			dLL_dV = data.Y - D_lambda; 			% NumObs x 1
		end

		%% Compute grad
		grad = compute_XprimeY2(dims, data.Xparts, dLL_dV)';		% NumParams x 1
		clear dLL_dV;
		if takeNegative; grad = -grad; end;
		varargout{1} = grad;
	end
	
	if nargout >= 3
		%% Compute d2LL_dV2
		if modelId == 1  % Logistic-Binomial
			d2LL_dV2 = -data.N .* p .* (1-p);				% NumObs x 1
		end
		if modelId == 2  % Poisson
			d2LL_dV2 = -lambda; % NumObs x 1
		end
		if modelId == 3  % Exponential
			d2LL_dV2 = -D_lambda; % NumObs x 1
		end
		
		hessian = compute_XprimeY_X2(dims, data.Xparts, d2LL_dV2);	% NumParams x NumParams
		clear d2LL_dV2;
		if takeNegative; hessian = -hessian; end;
		varargout{2} = hessian;
	end
end
